Pair condensation in a Finite Trapped Fermi Gas 
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We study signatures of the superfluid phase transition and the pseudogap phase in a trapped finite- 
size system of spin-1/2 fermions with infinite scattering length. Applying the auxiliary-field quantum 
Monte Carlo (AFMC) method in the canonical ensemble, we compute the energy-staggering pairing 
gap, heat capacity, and condensate fraction as a function of temperature. Our calculations reveal 
clear signatures of the superfluid phase transition in all three quantities, including a signature of 
the lambda peak in the heat capacity. Comparing the temperature dependence of these quantities 
shows that the pairing gap does not exhibit an obvious pseudogap effect in this system. 
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Cold atomic Fermi gases are an important paradigm in 
the study of strongly interacting quantum systems and 
of superfluidity. They provide clean and experimentally 
tunable strongly interacting systems that can be used 
as testbeds for the development of many-body methods, 
and their study can provide insight into the behavior 
of such diverse systems as neutron matter, quark mat- 
ter and high-temperature superconductors [1-5]. They 
have been widely investigated experimentally, and many 
of their properties, such as thermodynamic functions, 
density profiles and spectral functions, have been mea- 
sured [4, 6]. 

A key goal in the study of cold atomic Fermi gases 
is to accurately describe the superfluid phase transi- 
tion in the unitary gas and, in particular, to under- 
stand to what extent pairing persists in the normal phase 
above Tc. In cold atomic Fermi gases, Bardcen-Cooper- 
Schrieffcr (BCS) superfluidity is continuously connected 
to Bosc-Einstein condensation (BEC) through a strongly- 
interacting crossover regime. This connection is en- 
abled through s-wave scattering length a, which com- 
pletely characterizes a zero-range interaction. In the BCS 
regime, where a is small and negative, pairing and super- 
fluidity occur simultaneously at a critical temperature Tc- 
In the BEC regime, where a small and positive, pairing 
occurs prior to condensation at a temperature T* ^ Tc 
owing to the existence of a deep two-body bound state; 
in this regime a nonzero gap exists above Tc- In the 
crossover regime, however, the effects of pairing are not 
completely understood. 

Pairing is widely believed to occur above Tc in the 
crossover regime, and in particular in the unitary limit of 
infinite scattering length, giving rise to a so-called "pseu- 
dogap" phase (see Ref. [1] for an introduction). But the 
experimental characterization of this phase is incomplete. 
New radiofrequency spectroscopy techniques have pro- 
vided evidence for pairing at and possibly above Tc by 
observing back-bending in the dispersion relation [7-10] 
(there is some ambiguity for T > Tc because of the pres- 
ence of a universal tail [11] in the momentum distri- 
bution). However, experiments have not yet determined 
the complete temperature dependence of the pairing gap 



across the phase transition, which could provide one of 
the key signatures of the phase. Recently the heat capac- 
ity and compressibility were measured across the phase 
transition [12]; surprisingly, these quantities displayed no 
obvious signatures of a pseudogap above Tc, a result con- 
sistent with measurements of the equation of state for the 
pressure from another set of experiments [13-15]. 

Numerous Monte Carlo studies have been performed 
to study cold atomic Fermi gases. However, only a 
few [14, 17-19] have attempted to address the existence 
of a pseudogap phase above the superfluidity critical tem- 
perature. Even fewer have considered the thermodynam- 
ics of the trapped gas [20] or of finite-size systems. 

Here we study signatures of the superfluid phase tran- 
sition in the trapped, two-species ("spin- up" and "spin- 
down" ) , finite-size Fermi gas in the unitary limit of in- 
finite scattering length. We apply the auxiliary-field 
Monte Carlo (AFMC) method in the canonical ensem- 
ble (of fixed number of spin-up and spin-down atoms) to 
accurately calculate the pairing gap, condensate fraction, 
and heat capacity as a function of temperature. We de- 
fine the pairing gap as one-half the energy required to 
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l,Ni - l)]/2, where E{N-i,Ni) is the energy of 
the system with A'^ spin-up and A'^j, spin-down atoms [6] . 
This is a model-independent quantity. Our calculations 
provide the first ab initio results for the temperature de- 
pendence of these quantities in a finite, trapped atomic 
system, and the first ab initio calculations for the energy- 
staggering pairing gap and heat capacity across the su- 
perfiuid phase transition in any system of cold atoms. 
Our results are particularly relevant in light of recent ad- 
vances in experiments with finite-size quantum gases [21]. 

Hamiltonian and model space. — We consider equal 
numbers A'^ = = N/2 = 10 of two species of fermions 
interacting at very short range in a spherical harmonic 
trap with frequency w. The interaction V{r) is modeled 
as a contact interaction V{r) = Vo6{r), which acts only 
in the s-wave channel and therefore allows only parti- 
cles of different species to interact. We apply the AFMC 
method in the framework of a configuration-interaction 
(CI) shell model space spanned by the single-particle 
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eigenstates \nlm) of the harmonic trap (n is the radial 
quantum number, I is the orbital angular momentum and 
m is its projection) with energies = (2n + / + 3/2)huj. 
The model space is truncated to A/'max oscillator shells, 
i.e., we take only oscillator states with 2n + I < TVnax- 
We denote the total number of such states by TVs. By 
studying the convergence of observables as a function of 
TVmax we determine a cutoff for which the observables 
are well- converged at all temperatures of interest. The 
interaction strength Vb for each value of TV„iax is tuned 
to reproduce the exact ground state energy Eq of the 
two-particle system (Eq ~ 2huj at unitarity). 

AFMC method.— The AFMC method works by apply- 
ing a Hubbard-Stratonovich (HS) transformation [22] to 
the propagator e^*^^ (where /3 is the inverse temperature 
T) and evaluating the resulting multidimensional path in- 
tegral by Monte Carlo methods. We apply the CI shell 
model formalism commonly used in nuclear physics [23- 
26] . The thermal propagator is expressed as a path inte- 
gral 

e-^" « J D[a]G,U, , (1) 

where the integration is over auxiliary fields (7(t) that 
depend on the imaginary time t (0 < r < [/q- is 
the many-particle propagator for a system of nonintcr- 
acting fcrmions moving in external fields cr(T), and Go- 
is a Gaussian weight. To facilitate numerical evaluation 
the imaginary time is discretized into a finite number 
of points r„ = nA/?, where A/3 = ji/Nr- The calcula- 
tions become exact in the limit A/3 — 0. We find that 
observables scale linearly versus A/3 for A/3 < 1/32 and 
take the limit A/3 — > by a linear extrapolation. 

Observables. — In AFMC, a set of configurations ct^ are 
sampled according to the distribution Wa = Go-lTrt/crl 
using the Metropolis algorithm [25]. Thermal expecta- 
tion values of observables are then calculated from 

Tr(e-/5ff) ^ Efe*.. ' ^ ' 

where {0)cr = Tr(C'C/cr)/TrC/o- is the expectation of an 
operator O at a given configuration crj. of the auxiliary 
fields and = Tt{U„)/\TtU„\ is the Monte Carlo sign. 

Unlike most other AFMC calculations, we calculate 
the traces in (2) and in the weight function Wa at fixed 
particle numbers Nf,N^^, i.e., in the canonical ensemble. 
This is accomplished by using an exact particle-number 
projection via a discrete Fourier transform [25] . 

Monte Carlo sign. — In general the sign $o- is a sample- 
dependent complex phase. If the fluctuations in $0- be- 
come comparable to its average value, the statistical er- 
ror in (2) becomes excessively large, giving rise to the 
so-called Monte Carlo sign problem. However, attractive 
contact interactions are known to have good sign ($0 — 1 
for all a) in the grand-canonical ensemble [27], thereby 



avoiding the problem. This is also the case in our canon- 
ical ensemble calculations, and can be seen as follows. 
The sign is determined by particle-projected two-species 
trace 

Ttn^nA = {t^n.uI) {ttn^U^) , (3) 

where (U^) sltc the auxiliary-field propagators for 
spin-up (spin-down) particles. For an attractive contact 
interaction, and C/j: are time-reversal invariant, and 
both particle-projected traces on the right-hand side of 
(3) are real. For the spin-balanced system ^ N^, we 
have TrAT^^AT^C/cr = {Ttn^U^)'^ > so the canonical sign 
is also 1. 

Accuracy and convergence. — The accuracy of our cal- 
culations is determined by two factors. First, the contact 
interaction is not an exact representation of the usual reg- 
ularized delta function V{r) ~ 6{r)d/dr used to model 
cold atoms [28]. We checked its accuracy against the 
known analytic solution of the three-particle system and 
found the energy to be accurate within 1%. 

Second, for the calculations to be accurate, the model 
space must be large enough to account for both interac- 
tion effects and thermal excitations. We have accounted 
for this by using model spaces that are sufficiently large 
for the observables of interest to be well- converged for the 
relevant particle number and temperatures. We demon- 
strate the convergence in the insets to Figs. 1(a) and 1(b); 
the pairing gap Agap is well-converged by TVmax = 9, and 
the condensate fraction n by TVmax = 11. We used these 
values of TVmax in our calculations, along with A/'max = 11 
for the heat capacity (which is also well-converged). 

Algorithmic improvements. — At large /3, the matrix 
representing the one-body propagator becomes ill- 
conditioned and requires numerical stabilization. How- 
ever, the computational effort for the usual stabilization 
method [29] increases from 0{Mg) to 0{Af^) when gen- 
eralized from the grand-canonical ensemble to the canon- 
ical ensemble [30]. We have devised an improved algo- 
rithm for stabilization in the canonical ensemble which 
scales as 0{Afg) [31] and thereby makes the calculations 
tractable for the large model spaces necessary for conver- 
gence in this system. In particular TVmax = 11 contains 
TVs = 364 single-particle states and to our knowledge 
is the largest model space used to date for canonical- 
ensemble AFMC calculations. 

Pairing gap. — We calculated the pairing gap Agap 
by particle- number reprojection [32], in which only one 
Metropolis walk is needed to calculate all three energies 
E{N-^ - 1, A^; - 1), E{N^,Ni - 1) and E{N-^,Ni) con- 
tributing to Agap. We sampled the fields according to 
the distribution Ga\TrN.f,N^-iUa-\- 

Our result for the pairing gap is shown in Fig. 1(a) 
as a function of temperature. Here Agap is plotted in 
units of ep = 4.0 /iw [33] and the temperature is given 
in units of Tp = ep/kB, where fc^ is the Boltzmann 
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FIG. 1. Signatures of the superfluid phase transition in the 
trapped spin-balanced TV = 20 atom system (A'^ = A''^ — 10). 
(a) AFMC resuhs for the pairing gap Agap vs. temperature. 
The inset shows the convergence of the gap as a function 
of the number of oscillator shells A/lnax at T/Tf = 0.2. (b) 
Condensate fraction vs. temperature. Solid circles: the scaled 
occupation no of the T = pair wavefunction. Open squares: 
the largest scaled eigenvalue n of the L — pair correlation 
matrix (4). Solid line: noninteracting result. The dashed 
horizontal line indicates the noninteracting upper bound of 
l/{N/2). The inset shows the convergence of n as a function 
of A/'max at T/Tf — 0.125. (c) Heat capacity vs. temperature. 
The AFMC results (solid circles) are compared with the heat 
capacity of noninteracting fermions in the trap (solid line). 
The vertical dashed line in the three panels corresponds to a 
temperature of T/Tf = 0.175 (see text). 



constant. At high temperature, the finite size of the sys- 
tem causes Agap to remain slightly above zero. As the 
temperature decreases, the gap begins to depart from 
its high-temperature behavior at around T /Tp « 0.175; 
it then rises rapidly before saturating at approximately 
T/Tp = 0.07, providing a clear signature of a superfluid 
phase transition. We estimate the zero-temperature pair- 
ing gap by averaging the values at the lowest two tem- 
peratures, obtaining Agap = 0.271 (32)£i7. This result 



is consistent with fixed-node diffusion Monte Carlo [34] 
and density functional theory [35] calculations. 

It has been predicted that at zero temperature the 
energy-staggering pairing gap in a trap will be suppressed 
compared to its value in the homogeneous system, since 
in a trap unpaired particles can locate themselves in the 
edge of the cloud where the energy required to add an 
extra particle to the system is smaller [36]. More pre- 
cisely, Agap was predicted to scale at zero temperature 
as Agap oc (TV + If/^fujj. Using ep = {iNf/^Huj [6], 
one therefore expects Agap/eF oc TV~^/^ for large TV. 
We compare this scaling with our results by comput- 
ing the pairing gap at low temperature for TV 10 
(TVt = = 5). We find Agap(5,5) = 0.294(21), 
so that Agap(10, 10)/Agap(5,5) = 0.92(13), surprisingly 
close (for such small TV) to the expected value of 0.866, 
although results at larger values of TV would be required 
for a conclusive comparison. 

Condensate fraction. — We determined the condensate 
fraction from the two-body density matrix. For uniform 
systems, there is an equivalence between off-diagonal 
long-range order (ODLRO) and the existence of a large 
eigenvalue in the two-body density matrix [37]. Thus, 
the existence of a large eigenvalue in this matrix is con- 
nected with the presence of superfluidity. In the trapped 
system, the orbital angular momentum L is conserved, 
and we calculate the pair correlation matrix [38] 

CL{ab,cd) = {Alj^^^^{ab)ALMnicd)) > (4) 

where (•) denotes a thermal expectation value as in 
Eq. (2). Here A]^j^,j^^{ab) is a pair creation operator for 
two particles in orbitals a = (n^, la) and b = (ni,, 4) cou- 
pled to total angular momentum L with magnetic quan- 
tum number M, with the first particle being of the t 
species and the second the l species. The expectation 
values on the right-hand side of (4) are independent of 
M because of rotational invariance. Up to a (2L -t- 1)- 
fold degeneracy, the eigenvalues of Cl {ab, cd) are exactly 
those of the matrix {a\^^^al^^^adm.diOLcra,^) , which is 
the the orbital-space form of the density matrix studied 
in the uniform gas (via ODLRO) [27, 39-41]. We find 
that the largest eigenvalue always occurs for L = 0. 

As the temperature decreases, the maximum eigen- 
value Amax becomes much larger than all others 
and satisfies constraints which allow for two natu- 
ral definitions of a condensate fraction. The cor- 
responding eigenvector Lp{ab) defines a pair creation 
operator = X)a b 'i'*('^^)^oot4.('^^) '^hich satisfies 
(B^ B) = Aniax, indicating that Amax is the occupa- 
tion of a two-body wavefunction. For free fermions, 
< Amax £ 1, so no morc than two fermions can oc- 
cupy a paired state. However, in the presence of inter- 
actions Amax can exceed 1, effectively allowing fermion 
pairs to "condense" . One can then show that Amax sat- 
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isfies < Aniax < B{Afs,N), where 

B{Af„ N) = NiAfs - N/2 + 1)/2A(, < N/2 (5) 

and TVs is the number of single-particle orbitals for a sin- 
gle species [42]. In the limit of large Ms, 3(^/3, N) ap- 
proaches N/2. In fact for TVmax = 11 oscillator shells 
TVs = 364, and B(TVs,TV) = 9.75 is quite close to 
TV/2 = 10. Moreover, Amax = B{Afs,N) can be achieved 
for a particular many-body wavefunction (which is not 
necessarily an eigenstate of the system) . We may there- 
fore define a condensate fraction by n = A,„ax/(TV/2), 
where < n < 1. 

We can also define a condensate fraction in terms of 
the scaled occupation of the zero-temperature pair Bq = 
^■f(r = 0), i.e., no = (bIBq) / {N/2). In contrast to n, 
no measures the fraction of fully condensed pairs at finite 
temperature, ignoring contributions from paired particles 
not occupying the zero-temperature pair wavefunction. 
In general no{T) < n(T) and equality holds in the limit 
T = 0. 

Fig. 1(b) shows our result for ng and n as a function 
of temperature. The two condensate fractions are nearly 
identical for temperatures below T/Tp = 0.175, while 
TT-o < n at higher temperatures. For comparison, we plot 
n for a noninteracting gas (solid line) and its largest pos- 
sible value of 2/iV corresponding to {B^ B) = 1 (dashed 
line) . In both cases we see a rapid increase of the conden- 
sate fraction below T/Tp w 0.175, signifying the onset of 
a superfiuid state. 

The condition (B^Bo) > 1 provides a definite upper 
bound Tub in temperature for condensation to be present, 
since condensation requires the pair wavefunction to be 
occupied by more than one Fermion pair. In our system, 
this yields T^\,/Tf 0.175, a temperature similar to the 
transition temperature scales for the pairing gap and the 
heat capacity (see below for the latter). 

Pseudogap effects. — An effective way to assess pseudo- 
gap effects in this system is to compare the temperature 
dependence of the condensate fraction with the temper- 
ature dependence of the pairing gap. In the absence of 
other mitigating effects, preformed pairs of attractively 
interacting fermions would cause the gap to lead the con- 
densate fraction as the temperature decreases toward Tc. 
Comparing Figs. 1(a) and 1(b) indicates that this does 
not occur: First, Agap(Tiib) = 0.028(2), an insignificant 
fraction of its T = value. Second, below T/Tp « 0.175, 
the condensate fractions n and ng both grow simultane- 
ously with the pairing gap A. Thus, we conclude that in 
the TV = 20 finite-size system, Agap does not display any 
pseudogap effects. 

Similar to our results, the t-matrix calculations of 
Ref. [43] for a trapped gas found a pseudogap temper- 
ature only slightly higher than the critical temperature 
for the unitary trapped gas. In contrast, Monte Carlo 
calculations for a uniform gas found a nonzero gap A 



at temperatures significantly above T^ [17, 19]. In those 
calculations A was extracted from the spectral weight 
function (by fitting a BCS-like dispersion) rather than 
from energy staggering. 

Note that our conclusion in the finite-size system, 
where the phase transition is smoothed, does not imply a 
pseudogap effect in Agap will be absent for larger systems. 
At zero temperature, fixed-node quantum Monte Carlo 
and density-functional calculations predict that finite- 
size effects become negligible for more than 50 parti- 
cles [44], while lattice Monte Carlo computations have 
found that shell effects still persist at the 2% level above 
40 particles [45]. It would be interesting to use our 
method to study the finite-temperature behavior of the 
system as a function of TV. 

Heat capacity. — Fig. 1(c) shows our results for the 
heat capacity C = dE/dT as a function of tempera- 
ture. We calculated C by numerically differentiating 
the energy inside the path integral; this method takes 
into account correlated errors and greatly reduces the 
statistical error compared to differentiation after cal- 
culating E{T) [46]. At high temperatures, the heat 
capacity agrees with that of a noninteracting trapped 
Fermi gas (solid line). As the temperature decreases 
below T/Tp = 0.175, it begins to deviate from its 
non-interacting limit, remaining elevated until dropping 
rapidly at around T/Tp = 0.11 — a smoothed but ob- 
vious signature of the lambda peak observed recently in 
the thermodynamic limit [12]. The temperature at which 
this structure emerges, T/Tp « 0.175, is commensurate 
with the rise of the pairing gap and the temperature at 
which the condensate fraction n exceeds its noninteract- 
ing limit. 

In conclusion, we have calculated the energy stagger- 
ing pairing gap, condensate fraction, and heat capacity 
for an unpolarized system of TV = 20 trapped cold atoms 
at unitarity using a quantum Monte Carlo method. Clear 
signatures of the superfiuid phase transition exist in each 
of these quantities. In this finite-size trapped system, the 
gap does not lead the condensate fraction as temperature 
decreases and thus provides no clear signature of a pseu- 
dogap phase. The same method can be applied on either 
side of the BEC-BCS crossover, and could potentially be 
used to study finite-size scaling. 
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